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The nonlinear optical Galton board (NLOGB), a quantum walk like (but nonlinear) discrete 
time quantum automaton, is shown to admit a complex evolution leading to long time thermalized 
states. The continuous limit of the Galton Board is derived and shown to be a nonlinear Dirac 
equation (NLDE). The (Galerkin truncated) NLDE evolution is shown to thermalize toward states 
qualitatively similar to those of the NLOGB. The NLDE conserved quantities are derived and used 
to construct a stochastic differential equation converging to grand canonical distributions that are 
shown to reproduce the (micro canonical) NLDE thermalized statistics. Both the NLOGB and the 
Galerkin-truncated NLDE are thus demonstrated to exhibit spontaneous thermalization. 


I. INTRODUCTION 


At the fundamental level, quantum theory is linear. 
Yet, nonlinear models are often useful to take into ac¬ 
count interaction in an effective manner. Two exam¬ 
ples are the so-called nonlinear optical Galton board 
(NLOGB) [I] and the nonlinear wave equations de¬ 
scribing the dynamics of Bose-Einstein condensates 
(BEC). Though the NLOGB is discrete and wave 
equations are by definition continuous, these models 
have much in common. Indeed, the NLOGB is essen¬ 
tially a nonlinear quantum walk (QW), and the for¬ 
mal continuous limits of linear QWs are wave equa¬ 
tions [Ml- Typical such wave equations are the Dirac 
or the Schrodinger equation which nonlinear version, 
called the Gross-Pitaevskii equation (GPE), is used to 
model BEC [3- Also, QW descriptions of BEC have 
been proposed in HUD- 

Finally, numerical solutions of continuous wave 
equations are actually solutions of discrete systems ap¬ 
proximating the continuous equations. 

The NLOGB can be seen as a discrete model of non¬ 
linear waves similar to those which propagate in BEC. 
One can therefore expect the NLOGB to display prop¬ 
erties similar to those of the standard nonlinear model 
of BEC: the GPE. One such property which has until 
now never been explored on the NLOGB nor, more 
generally, in the context of QWs and quantum au¬ 
tomata, is the so-called spontaneous thermalization. 

In the context of (nonlinear) BEC, microcanoni- 
cal equilibrium states are well-known to result from 
long-time integration of the so-called truncated (or 
Galerkin-projected) Gross-Pitaevskii equation (GPE) 
and involve a condensation mechanism [Z1II0HI2]- Fur¬ 
thermore, such thermalization is also known to happen 
in discretized (rather than spectrally-truncated) GPE 
[151 . Classical Galerkin-truncated systems have been 


studied since the early 50’s in fluid mechanics. In this 
context, the (time reversible) Euler equation describ¬ 
ing spatially-periodic classical ideal fluids is known 
to admit, when spectrally truncated at wavenumber 
fcmaxj absolute equilibrium solutions with Gaussian 
statistics and equipartition of kinetic energy among 
all Fourier modes nunj. Furthermore, the dynam¬ 
ics of convergence toward equilibrium involves a direct 
energy cascade toward small-scales HHHS. 

The aim of the present work is to study thermal¬ 
ization phenomena in a spatially-periodic version of 
the NLOGB and relate it to the thermalization of its 
(Galerkin-truncated) continuous limit. 

The paper is organized as follows. Section [D] is de¬ 
voted to the definition of the NLOGB model and its 
numerical solution. The main result of this section 
is to display and characterize the complex behavior 
of the log-time regime. Section [Hi] is devot ed to the 
behavior of the continuous limit (sect. Ill A I, its con¬ 
served quantities (sect III B) and the long-time behav¬ 
ior and thermalization of its Galerkin-truncated ver¬ 
sion (sect. |III C ). Finally section IV is our conclusion. 
Technical details are given in appendices. 


II. NONLINEAR DISCRETE TIME 
QUANTUM WALK 

A. Fundamentals 


Consider a quantum particle endowed with an in¬ 
ternal degree of freedom and a lattice on which this 
particle can move in discrete time. A Discrete Time 
Quantum Walk (DTQW) is an automaton which con¬ 
ditions the motion of the particle on the state of its 
internal degree of freedom [2(3]. Let us remark that 
whilst a continuous-time version of QW (CTQW) - 
living in continuous time and discrete space - has been 
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introduced in the literature [21] , we will not deal with 
it in the present work. In this article, we focus on a 
discrete time nonlinear quantum walks (DTQW) de¬ 
fined on the discrete circle and on particles described 
by a two components complex wave function. The dis¬ 
crete time spatially periodic quantum walk is defined 
by the following equations: 
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The index m = 0,...,N — 1 labels points on the dis¬ 
crete circle and the index j £ N labels instants. At 
each time j and each point m, ipj m are the two com¬ 
ponents of the wave-function on a certain space- 
and time-independent spin basis (&_ The above 

finite difference equations were derived in [I], albeit 
for QWs on the unrestricted line m The parameter 
g fixes the strength of the nonlinearity. For g = 0, 
equations ([!]) coincide with the evolution equations of 
the standard Hadamard walk. The particle number at 
time j 


n i = E (i^-j 2 + iv^J 2 ) = E( 2 ) 

m m 

is independent of j i.e. it is conserved by the walk and 
normalized to 1. We will henceforth denote it by II. 


B. Asymptotic behavior of the DTQWs 
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As displayed in FigJT^i, the family of DTQWs defined 
by equations ([l]) exhibits a very complex dynamics, 
much richer than the dynamics of the Hadamard walk 
shown for comparison in Figjljx 

Of particular interest is the j —> oo asymptotic 
statistics of the particle number spatial distribution 
II j,m = \*p-J 2 + |^y m | 2 - Let A p be a positive real 
number much lower than unity and compute, at all 
times j, the proportion Hj(p)Ap of values of the posi¬ 
tion m for which Hj m lies in (p,p+ Ap). The quantity 
Hj ( p ) is a discrete equivalent of the probability density 
function (PDF) of Hj m at time j. 

Direct numerical simulations (DNS) (see Fig[l]a) al¬ 
low to directly estimate Hj ( p ). It is found that Hj ( p ) 
tends towards a stationary distribution H*{p) which 
depends on the initial condition. Figure Fig. (|2ja) 
displays how the particle number Hj, m typically de¬ 
pends on m at fixed large values of j and Fig. @b) 
displays H*{p), estimated by computing histograms as 
explained above. 

The existence of H* ( p ) is typical of nonlinear 
chaotic systems. These systems also exhibit a 
great sensitivity towards initial conditions, and this 
sensitivity is confirmed by DNS of the NLOGB. 
Indeed, starting a DNS of the Hadamard walk with 


Figure 1. (Color online) Density plot of the time evolu¬ 
tion of the particle number density IIy m , as defined in 
equation for (a) the Hadamard DTQW (g= 0) and 
the NLOGB (<?=10 7r), with a symmetric initial condition 
+ ib+). Number of grid points N = 512. 


a symmetric initial condition delivers a numerical 
solution which is symmetrical at all times, whereas 
using the same initial condition in a DNS of the 
NLOGB delivers a numerical solution which is not 
symmetric (see Fig. [lja). This symmetry breaking 
becomes greater with the time j (see Fig. @b) 
and depends on the resolution of the DNS and the 
strength of nonlinearities. In particular Fig|3]b shows 
that the symmetry breaking starts from the round-off 
noise m that is of order 10 15 in our simulations. 
We have checked (data not shown) that adding to the 
initial condition a non-symmetric noise larger than 
the round-off noise produces the same growth rates 
for the symmetry breaking, but starting at the higher 
level of the added non symmetric noise. This confirms 
that the symmetry breaking is due entirely to the 
round-off noise. 
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in Appendix [A] The formal continuous limit of the 
NLQWs read: 
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Figure 2. (Color online) (a) Log-Lin plot of the particle 
number density IIj, m (see eq. @) at time j = 5x 10 5 (yellow 
cercles) for the NLOGB with the same initial condition as 
Fig. [l](black point) and for g= 10 n. Number of grid points 
N = 64. (b) Log-Lin PDF H(p) (see text 2nd paragraph 
in Sec. B) of ^lj,m at time j = 5 x 10 s for the NLOGB with 
the same initial condition as Fig. [l]and g= 10 7r. Number 
of grid points N = 64. 


III. NONLINEAR DIRAC EQUATION 

A. A nonlinear Dirac equation as continuous 
limit of the DTQWs 

The asymptotic aspects of the NLOGB dynamics 
can be understood by investigating the continuous 
limit of these walks. The method employed is the same 
as in [3 HOES and detailed computations are given 
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where I is the identity, 
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are the second and third Pauli matrices. The con¬ 
tinuous limit of the NLOGB is thus described by a 
nonlinear Dirac equation (NLDE). The nonlinearity is 
confined to the mass term, which depends quadrat- 
ically on the spinor T. Note that (spatially Two- 
dimensional) NLDE have also been used to describe 
experimental BEC on 2D hexagonal lattice pirl - ESl . 

The NLDE 0 is formally equivalent to Nambu- 
Jona-Lasinio-like equations (NJLE) (Nambu and 
Jona-Lasinio, P5h33| ) in 1+1 dimension, which de¬ 
scribe a nonlinear interaction between fermions with 
chiral symmetry. The constant g corresponds to a non¬ 
linear coupling constant and if g = 0, © degenerates 
into the Weyl equation. 


As detailed in Appendix [B] the validity of the 
continuous limit is best confirmed by using Fourier 
pseudo-spectral methods ED, which are precise and 
rather easy to implement. In particular, Fig. 0 dis¬ 
plays for different values of g the relative difference 
between the solution of equations ([I]) and 0 as a func¬ 
tion of the e parameter which controls the continuous 
limit. 


Fig. © shows the typical profile of the asymptotic 
n(T, X) and the stationary distribution H(p ) of this 
density, as obtained from a Galerkin-truncated sim¬ 
ulation of the NLDE, de-aliased in a way that en¬ 
sures conservation laws in the truncated system, see 
Appendix [0 Both plots are strikingly similar to the 
corresponding plots presented in Fig. [2]a and 0b 
obtained by numerically integrating the NLOGB. In 
other words, the NLOGB and the Galerkin-truncated 
NLDE seem to have very similar asymptotic behav¬ 
ior. We will now analyze in detail the asymptotic be¬ 
havior of the Galerkin-truncated NLDE. We will first 
identify the conserved currents for the NLDE (Section 


III BI and then show that the asymptotic statistics 


Galerkin-truncated NLDE is identical to the so-called 
grand canonical statistics (III C). 
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Figure 3. (Color online) (a) Density plot of the time evo¬ 
lution of the particle number density Tlj,m (see eq. <§) 
with the same conditions as Fig{T]b , but in the short 
time regime, (b) Asymmetry measure F j versus time, 
for different values of the number of grid points. L, 
= Z)m=o _1 “ E™=jv /2 n j.™ and same conditions as 

Fig^b. 


B. Lagrangian formulation and conserved 
quantities 

The NLDE derives from the following Lagrangian 
density: 

£(T, tft) = 1 ) - (d^)Y*] - | (M$) 2 , 

(7) 

with 



Figure 4. (Color online) Log-Log plot of the relative differ¬ 
ence Sj at time j = 400, defined as 

where Hqw is the particle number density If Jjm obey¬ 
ing to the NLOGB finite difference equations and II d 
is the particle number density TT(T), X rn ) obeying to the 
NLDE. The relative difference is shown for several values 
of e=27 v/N, (from right to left) N=2 6 ,2 7 ,2 8 ,2 9 . The initial 
condition is a symmetric Gaussian distribution (0, X m ) = 
^\/5 (b- + ib+), where the Gaussian shape /( X m ) = 
exp [—X^l\/2 <r 2 ) and a = 10 A*. 


N -1 ( n gn'- n fl)‘ 

m =° (n dWZT 
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^, 7 1 = i 02, 7 5 = i 7°7 1 , ’F = 
and do = dr : d\ = dx- 

There are two conserved currents and these gener¬ 
ate three integrals of motion (conserved quantities). 
The first current is simply the 2-current J^ = TyC ip 
associated to the U( 1) invariance of the NLDE. The 
corresponding integral of motion is the total particle 
number: 

II[’F,\E' t ] =Ju(T,X)dX, (9) 

where II (T,X) = 4/1 (T, A')T(T, X), which is usually 
normalized to 1 . 

The other current is associated to the space-time 
translation invariance of the NLDE and is the stress- 
energy tensor 

T'“ / ('F,’F t ) = -[$7^’F) - {d v ^)^^}-if v C, 

2 ( 10 ) 

where 77 ^" = diag(l,—1). The associated conserved 
quantities are the energy E and the momentum P, 
which are defined by 


7 = = 



N = 7o + ^75, 


( 8 ) £[)F,Tt] = J T 00 ('F(T,X),T t (T,X)) dX ( 11 ) 
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Figure 5. (Color online) (a) Log-Lin plot of n(T, X) at 
time T = 5 x 10 5 (red square) obeying the NLDE with the 
same initial condition of Fig|4] (black point) for g =107r. 
Number of grid points N = 64. (b) PDF H(p) of II. 


and 


C. Thermalization in the Galerkin-truncated 
NLDE 


If one studies the NLDE on the circle, it is natural 
to write at all times the spinor ^(T, X) as a spatial 
Fourier series and to replace the NLDE by an evolution 
equation obeyed by the time-dependent Fourier coef¬ 
ficients ^(T, k). In performing a Galerkin truncation 
13- r )| . one retains only a finite number of these coef¬ 
ficients as dynamical variables, say 'F(T', k) with k = 
— y,..., — 1, and replaces the exact NLDE dynamics 

by a new dynamics which, at small k, approximates 
at least formally the original NLDE dynamics. By 
Fourier transforming the ^(T, k), k = —jf,If — 1, 
back to original physical space (i.e. the circle), one ob¬ 
tains a set of N spinors \P m (T), m = 0..., N — 1, which 
are to be interpreted as the values ’L(T, X m ) taken 
by the spinor field ^(T, X) at point X m = (see 
Appendix |B| . The spinors ^(T, X m ) are on the same 
footing as the k = y — 1, and can 

be viewed as the dynamical variables of the Galerkin- 
truncated NLDE. We now denote by 'F(T) the collec¬ 
tion {'F m (T) = \I/(T, X m ),m = 0, ...IV — 1}. 

All integrals over space of quantities involving the 
Dirac field can be replaced by Riemann sums. Thus, 
the total particle number, the energy and the momen¬ 
tum can now be viewed as functions of the collection 
(’i r (T),’i'*(T)). These functions will still be denoted 
by II, P and E and are conserves by the Galerkin 
truncated dynamics, see Appendix [B| 

We now introduce two Lagrange multipliers /i and 0, 
define Hg^ = E — OP — pH and consider the following 
stochastic differential equations: 


d^f? m 

hr" 

dT 


dHg M 
9*rn 
dHep 


VCm(T) 


(15) 


where r) is a real coefficient and the £ m (T)’s are com¬ 
plex independent Gaussian white noises m with cor¬ 
relation functions 


< U(T)C(T') >= S mn S(T - T'). (16) 


P['f,'L t ] = fr»(*(T,n*HT,X))dX, (12) 


The density / of this stochastic process obeys the ex¬ 
act Fokker-Planck equation [ID EH EE] 


with 


T 00 (T,d' t ) 
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The stationary solution f.^g^ of this equation is the 
Gibbs state 


T 01 (^,#) = - l - [^7°(9a-^) - (9v v k)7° v k] • (14) 


■W = v — exp 



(18) 
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Figure 6. (Color online) PDF H(p) of the thermalized 
state density II(A') for the NLDE and for the stochastic 
equations (151 (GE) for g = 10-7T. The conserved quantities 
and the noise coefficient are: a) E= -19.8, P = 0, II = 1 
and r/ =1; b) E= -17.12, P = 0, II = 1 and r/ =1.6, at time 
r=312. 


which is the so-called grand canonical distribution 
with inverse temperature 2 /rj 2 . 

We have simulated the stochastic equations (151 to 
numerically obtain fields with statistical distributions 
corresponding to the distribution f v g M and compare 
it with the asymptotic long-time distribution of the 
Galerkin-truncated NLDE. 

The comparison between the thermalized fields cor¬ 
responding to the nonlinear Dirac equation and the 
Gibbs states corresponding to Eq.( 15) are displayed in 
Fig. [6]and confirm that the Galerkin-truncated NLDE 
and the stochastic equations © are described by very 
similar distributions (see Appendix |B|. 

Fig. [| has been generated in the following man¬ 
ner. We have first used © with rj = 1 and = 1.6 


to produce two Gibbs states with particle number II 
fixed to unity and vanishing momentum (see the end 
of appendix H- The energies of these two states are 
respectively E = —19.8 and —17.12. We then have 
generated initial data for the nonlinear Dirac equation 
([3]) with the same particle number (equal to unity) 
and energies. This has been done by multiplying the 
Gaussian initial field by a suitably chosen spatially de¬ 
pendent phase. 


IV. DISCUSSION 

We have considered the NLOGB confined to the cir¬ 
cle and we have shown that the continuous limit of 
this NLOGB is a NLDE identical to the NJLE-model. 
Pseudo-spectral numerical simulations reveal that the 
asymptotic behavior of the NLOGB is similar to the 
asymptotic behavior of the Galerkin-truncated NLDE 
and we have shown that the associated asymptotic 
statistics is identical to the grand-canonical statistics. 
Thus, both the NLOGB and the Galerkin-truncated 
NLDE exhibit spontaneous thermalization. Strictly 
speaking, the NLOGB on an unrestricted line, as pre¬ 
sented in |T], cannot be studied with the same method. 
In fact the comparison between the NLOGB and the 
Galerkin-truncated NLDE is possible only when the 
NLOGB admits a finite number of wave numbers. 

Previous work on other nonlinear quantum walk [ 35] 
suggests that this observed spontaneous asymptotic 
thermalization is not a particular feature of the sys¬ 
tems studied in this article, but will also be encoun¬ 
tered in other nonlinear quantum walks, whatever the 
dimensions of the underlying physical space or of the 
coin space may be. It is obvious that quantum walks 
which thermalize will explore space in a very differ¬ 
ent manner from walks which do not thermalize, and 
their importance for quantum computing should cer¬ 
tainly be explored in depth. In a different direction, 
it would be interesting to exhibit and analyze sponta¬ 
neous thermalization in QWs couple to synthetic gauge 
fields pi [251 HD] , 


Appendix A: Derivation of continuous limit 

Consider for all (re, j) £ N 2 , the collection W” = 
(’P k,m)k=nj,mei,- This collection represents the state 
of the NLOGB at ’time’ k = nj. For any given re, the 
collection S n = (WC’OjeN thus represents the entire his¬ 
tory of the NLOGB observed through a stroboscope of 
’period’ re. The evolution equations for S n are those 
linking W™ +1 to WJ 1 for all j. The method employed 
here to obtain the continuous limit of a generic S n was 
introduced in PI 12S]. 

One first introduces a time-scale r, a length-scale A, 
an infinitesimal e and interpret the space-index rre as 
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referring to position x m = me A = mAx and the time 
index j as referring to the instant tj = jer = jAt. 
The formal continuous limit is obtained expanding the 
equations defining S n in Taylor series around e = 0 
and by letting e tend to zero. For the limit to exist, 
all zeroth order terms of the Taylor expansion must 
identically cancel each other and the differential equa¬ 
tion describing the limit is then obtained by equating 
to zero the non identically vanishing, lowest order con¬ 
tribution. 

The original NLOGB S 1 does not admit a continu¬ 
ous limit because the zeroth order terms do not can¬ 
cel each other identically. The equations defining S 2 
read: 

tp (tj + '2At, x m ) — — [tF[<p (tj , x m + As)] + (tj , x m Ax')] 

ip^ (tj + 2 At, x m ) — — [J-[(p (tj , Xm — Ax)] — tF[cp+ (tj , Xm + Ax)] 

where 

= e‘sr(t j ,»m+A»)|^-( ( . iIm + Al ) i . 

e ig|l/.-(t,,-,* m +Ax)| 2 jp+( tjjXm + Ax) (Al) 

and 

FmjjXm)} = e ia ^' x ^<t>{tj,x m ). (A2) 


These equations admit a formal continuous limit, 
which reads: 

(Idr-Vdx - 3 |'Ad('F,T t ))'F =0 (A3) 

where 

Ad (T, 'F 1 ') = 'F 1 Al'F, (A4) 

M = I~Y (A5) 

and T = t/r and X = x/X. 

The operator P is self-adjoint and its eigenvalues 
are —1 and +1. Two eigenvectors associated to these 
eigenvalues are 

cos^) 6_ + (sin0 6+ (A6) 




and 



The family ( B _, B+) forms an orthonormal basis of the 
two dimensional spin Hilbert space. In this new basis, 


equation (A3) reads: 


(l&r-vadx- ^-M(>F,'F t ))'F = 0 (A8) 


where 


Appendix B: Numerical Methods 


We restrict ourself to 27r-periodic boundary con¬ 
ditions. A generic field ’F(A') is thus evaluated 
on the N collocation points X m = 2-irm/N, with 
m = 0,N — 1 as = f(I m ). The dis¬ 

crete Fourier transforms are standardly defined as 
T(A m ) = YXk=-N/ 2 ex P(ikXm)' 4 >k and the inverse 

i>k = jf J2mZ o exp (- ikX m ). These sums can 

be evaluated in only Nlog(N) operations by using 
Fast Fourier Transforms (FFTs). Spatial derivatives 
of fields are evaluated in spectral space by multi¬ 
plying by ik and products are evaluated in physi¬ 
cal space. The original QW equations can also be 
simply cast in this setting, as the translation opera¬ 
tor 'F m IWi is represented in Fourier space by 
Tfc —> d'fc exp (±ik2n/N). In this setting, the continu¬ 
ous limit is automatically taken when N is increased. 
As we can observe in Fig. Q the relative difference 
scales as expected as e for different values of ui. 
However the pseudo-spectral code solving the 
NLPDEs generates a problem called aliasing 52 , 
which means that high &-modes alias the amplitudes 
at lower fc-modes of the field. In that case the DFT 
is aliased and in general the fields needs to be de- 
aliased by proper spectral truncation. Here, we used 
the so-called 2/3-rule in all our numerical schemes in 
the same way as done in reference 52 - De-aliasing 
is fundamentally important to preserve the conserva¬ 
tion of the Galerkin truncated nonlinear dynamics as 
we can observe in Fig. ([7]). Indeed, although it is 
straightforward to show that Eq.([3]) can be written 

dA = (B1) 
u ^ m 

(B2) 


and thus formally conserves the energy, it can be 
shown that exact conservation requires proper de¬ 
aliasing (see appendix of ref. (IT]). 

As displayed in Fig |6] the statistical distributions 
generated by the NLDE dynamics Eq.([3]) and by the 
stochastic equations © are really close and this can 
be justified on very general grounds. 

First, by construction, the stochastic equations (151 
generate the grand canonical distribution ( [18] ) that is 
controlled by the inverse temperature 2/rf and the 
Lagrange multipliers /r and 9. On the other hand, as 
the spectrally-truncated dynamics (|3| conserves n, P 
and E, its long time behavior should be described by 
the so-called micro canonical distribution 


A/f('F,'F t ) = 


(A9) f ~6{E- E in )S{U - U in )S{P - P in ). (B3) 


M = ,+ f 


(A10) 


that is determined by the values (.E in ,n in ,P in ) of the 
conserved quantities given by the initial condition 'Flu. 
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As is well-known (42] , under very general circum¬ 
stances both grand canonical and micro canonical dis¬ 
tribution yield similar statistical results (provided that 
the 2/? 7 2 and the Lagrange multipliers /x and 0 have 
values that correspond to Ei n ,H in ,Pi n ). Note that the 
effect of a fixed value of /x in Eq. (151 amounts, at 
each time-step, to an overall multiplication of the field 


by (1 + fj,dt). Thus if we want the final result to have 
a fixed value of the total particle number II this can 
be obtained setting /x to zero and, instead, renormal¬ 
izing the field to the desired value of particle number 
at each tilme step. FigjG] indicates that, in this case 
(and zero values for P and 6 ), both distributions yield 
identical results for density fluctuations. 
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